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ABSTRACT 

We describe Monte Carlo models for the dynamical evolution of the nearby globular 
cluster M4. The code includes treatments of two-body relaxation, three- and four-body 
interactions involving primordial binaries and those formed dynamically, the Galactic 
tide, and the internal evolution of both single and binary stars. We arrive at a set of 
initial parameters for the cluster which, after 12Gyr of evolution, gives a model with 
a satisfactory match to the surface brightness profile, the velocity dispersion profile, 
and the luminosity function in two fields. We describe in particular the evolution 
of the core, and find that M4 (which has a classic King profile) is actually a post- 
collapse cluster, its core radius being sustained by binary burning. We also consider the 
distribution of its binaries, including those which would be observed as photometric 
binaries and as radial-velocity binaries. We also consider the populations of white 
dwarfs, neutron stars, black holes and blue stragglers, though not all channels for blue 
straggler formation are represented yet in our simulations. 
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1 INTRODUCTION 

The present paper opens up a new road in the study of 
the dynamical evolution of globu lar clusters. We adopt th e 
■ Monte Carlo method of Giersz (|Giers j Il998l . l200ll . |2006t) . 
which in recent years has been enhanced to deal quite real- 
istically with the stellar evolution of single and binary stars, 
to study the dynamical history of the nearby globular clus- 
ter M4. An earlier version of the code had already been used 
to st udy the dynamical history of to Cen (|Giersz fc Heggie] 
2003), but at that time the treatment of stellar evolution 
was primitive and there were no binaries. The new code has 
been thoroughly tested on smaller systems, by comparison 
with A-body sim ulations and observatio ns of the old open 
star cluster M67 l|Giersz fc He ggic 2008). There we showed 
that the Monte Carlo code could produce data of a similar 
level of detail and realism as the best A-body codes. Now for 
the first time we consider much richer systems, with about 
half a million stars initially, which are at present beyond the 
reach of A-body methods. 

This paper has a place within a long tradition of the 
modelling of globular star clusters, but the place is a dis- 
tinctive one. First, we are not concerned with a static model 
of a star cluster at the present day, like a King model. We 
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are concerned with issues where the dynamical history of 
the star cluster is important, where static models are un- 
informative. Secondly, our aim is to construct a model of a 
specific star cluster, rather than trying to understand the 
general properties of the evolution of a population of star 
clusters. T his has been done befor e, and a brief history is 
outlined in iGiersz fc Heggie] (|2008l ). but the present work 
takes these efforts onto a new level of realism, in terms of 
the description of stellar evolution, and dynamical interac- 
tions involving binary stars. 

This problem is not easy. Not only is it necessary to 
use an elaborate technique for simulating the relevant as- 
trophysical processes, but it is necessary also to search for 
initial conditions which, after about 12Gyr of evolution, lead 
to an object resembling a given star cluster. By "resembling" 
we do not simply mean matching the overall mass, radius 
and binary fraction of a cluster, for example, for two reasons: 

(i) We have found that values found for data in the lit- 
erature are highly uncertain, and different sources are con- 
tradictory. These data are usually derived, in some model- 
dependent way, from such data as surface-brightness profiles 
and velocity dispersion profiles, and we prefer to compare 
our models directly with this data, and not with inferred 
global parameters. 

(ii) We have found that, even if one achieves a satisfac- 
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tory fit to these profiles, the model may give a very poor 
comparison with the luminosity function. 

From these considerations we conclude that a model which 
aims to fit only the mass and radius of a star cluster (say) 
may be very far from the truth. 

Tackling this difficult problem is not just interesting, 
however. We have been motivated by a number of press- 
ing astrophysical problems. For example, the two nearby 
star clusters M4 (the subject of this study) and NGC 6397 
(which we shall consider in our next paper in this series) 
have rather similar mass and radius, and yet one has a clas- 
sic King profile, while the other i s a well-studied exa mple of a 
cluster with a "collapsed core" ()Trager et al.l[l995h . Among 
possible explanations one may consider differences in the 
population of binaries, which are known to affect core prop- 
erties, or in tidal effects. Indeed the present paper will show 
that these two clusters may be much more similar than one 
would suppose from the surface brightness profiles alone. 

A second motivation for our work is our involvement in 
observational programmes aimed at characterising the bi- 
nary populations in globular clusters. What differences (e.g. 
in the distributions of periods and abundances) should one 
expect to find between the core and the halo? These issues 
are important in the planning of observations, and in their 
interpretation. 

The cluster M4 is the focus of much of this effort be- 
cause it is nearby, making it a relatively easy target for 
deep observational study. It was the firs t globular cluster t o 
yield a deep sequence of white dwarfs (Ric her et al.lll995l ). 
More recently it has been subjected to an i ntensive ob- 
servational programme by the P adova group (|Bedin et al.l 
l200ll . 120031 ; Undersoil et alj|2006l ). which includes searches 
for radial-velocity bin aries in the upper main sequence 
l|Sommariva et all2008l ). It also turns out to be a cluster 
which (we conclude) started with only about half a million 
stars, which facilitates the modelling. Along with the open 
cluster M67, M4 was chosen by the international MODEST 
consortium, at a meeting in Hamilton in 2005, as the focus 
for joint effort by theorists and observers, to cast light on 
its binary population and dynami cal properties. M67 has 
been modelled very successfully bv lHurlev et al.l (|2005f ). us- 
ing iV-body techniques, and this paper represent the first 
theoretical step in a similar study of M4. 

The paper is organised as follows. First, we summarise 
features of the code and the models, the data we used, and 
our approach to the problem of finding initial conditions 
for M4. Then we present data for our best models: surface 
brightness and velocity dispersion profiles, luminosity func- 
tions, the properties of the binary population, white dwarfs 
and other degenerate remnants, and the inferred dynami- 
cal state of the cluster. The final section summarises our 
findings and discusses them in the context of work on other 
clusters, including objects to which we will turn in future 
papers. 



2 METHODS 

2.1 The Monte Carlo Code 

The details of our simulation method have been amply de- 
scribed in previous papers in this series. Each star in a spher- 



ical star cluster is represented by its mass, energy and an- 
gular momentum, and its stellar evolutionary state may be 
computed at any time using synthetic formulae for single 
and binary evolution. It may be a binary or a special kind 
of single star that has been created in a collision or merger 
event. 

Neighbouring stars interact with each other in accor- 
dance (in a statistical sense) with the theory of two-body 
relaxation. If one or both of the participants is a binary, the 
probability of an encounter affecting the internal dynamics 
is calculated according to analytic cross sections, which also 
determine the outcome. This is one of the main shortcom- 
ings of the code, as these cross sections are not well known in 
the case of unequal masses, and also the possibility of stellar 
collisions during long-lived temporary capture is excluded. 

A star or binary may escape if its energy exceeds a 
certain value, which we choose to be lower than the energy 
at the nominal tidal radius, in order to improve the scaling o f 
the lifetime with N, as explained in lGiersz fc Heggig {2008). 
This is the second main shortcoming of the models, as it 
leads to a cutoff radius of the model that is smaller than the 
true tidal radius, and this lowers the surface density profile 
in the outer parts of the system. 

A difficulty in applying the Monte Carlo code to M4 is 
that it employs a static ti de, whereas the orb it of M4 ap- 
pears to be very elliptical l|Dinescu et al.|[l999h . We have to 
assume that a cluster can be placed in a steady tide of such a 
strength that the cluster loses mass at the same average rate 
as it would on its true orbit. Some support for this proce- 
dure c omes from iV-body modelling. iBaumgardt fc Makind 
(|2003h show that clusters on an elliptical orbit between 
about 2.8 and 8.5kpc dissolve on a time scale intermedi- 
ate between that for circular orbits at these two radii, and 
that the dissolution ti me scales in almost the same way with 
the size of the system. IWilkinson et al.l l|2003h show that the 
core radius of a cluster on an elliptic orbit evolves in very 
nearly the same way as in a cluster with a circular orbit at 
the time-averaged galactocentric distance. 

All other free parameters of the code (e.g. the coefficient 
of iV in the Coulomb logarithm) take the optimal values 
found in the above study. 



2.2 Initial Conditions 

The initial models are as specified in Table Q] Many of these 
features (e.g. the properties of the binaries, except for their 
overall abundance) were inherited from our modelling of the 
open cluster M67, and tho s e were in turn mainly drawn from 
the work of iHurlev et al.l (|2005l ) . Some of the parameters 
were taken to be freely adjustable, and this freedom was 
exploited in the search for an acceptable fit to the current 
observational data. 



2.3 Observational data and its computation 

Our first task was to iterate on the initial parameters of 
our model in order to produce a satisfactory fit to a range of 
observational data at an age of 12Gyr. The data we adopted 
are as follows: 

(i) Surface brightness profile: here we used the compila- 
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Table 1. Initial parameters for M4 



Table 2. Properties of M4 



Fixed parameters 



Structure 


Plummer model 


Initial mass function 1 


KrouDa, Tout. & Gilmore (1993) in the 




range [O.1,5O]M 


Binary mass distribution 


KrouDa et al. (1991) 


Binary mass ratio 


Uniform (with component 




masses restricted as for 




single stars) 


Binary semi-major axis 


Uniform in log, 2(Ri + R 2 ) 




to 50AU 


Binary eccentricity 


Thermal, with eigenevolution 




(Kroupa 1995) 


Metallicity Z 


0.002 


Age 


12Gvr (Hansen et al. 2004) . 


Free parameters 


Mass 


M 


Tidal radius 


n 


Half-mass radius 


rh 


Binary fraction 


fb 


Slope of the lower mass 




function 


a (Kroupa = 1.3) 



tion bv lTraeer et~aP (|l995T ). where the surface brightness is 
expressed in V magnitudes per square arcsec. 

(ii) Radial veloci ty profile: this came from 
l|Peterson et al.l Il995l ). and is the result of binning 
data the on the radial velocities of nearly 200 stars. Strictly 
we should refer to the line-of-sight velocity dispersion, as 
"radial" velocity has a different meaning in the Monte 
Carlo model. 

(iii) The V- luminosity function (|Richer et al.|[20o3 . from 
which we considered the results for the innermost and outer- 
most of their four annuli). These data are lack correction for 
completeness, though th e completenes s facto r for the outer- 
most field is plotted in iRicher et al.l {2002). For main se- 
quence stars it is almost 100% down to V = 15, and drops 
steadily to less than 50% at V = 17. For the innermost field 
the completeness correction would be larger. 

Now we consider how to compare this data with the 
output of the Monte Carlo code. This includes a list of each 
particle in the simulation, along with data on its radius, ra- 
dial and transverse velocities and absolute magnitude, and 
numerous other quantities. To construct the surface bright- 
ness, we think of each particle as representing a luminous 
shell of the same radius, and superpose the surface bright- 
ness of all shells. While this procedure involves a minimum 
of effects from binning, or randomly assigning the full po- 
sition of the star, a shell has an infinite surface brightness 
at its projected edge. The effect of this, especially from the 
brightest stars, will be seen in some of the profiles to be 
presented in this paper. A similar problem arises in actual 
observations, and is often handled by simply removing the 
brightest stars from the surface brightness data presented. 
The output from the model is corrected for extinction (Table 
[3 The distance of the line of sight from the cluster centre 
is converted between pc (as in the model) and arcsec (as in 
the observations) using the distance in the same table. 



Distance from sun 
Distance from GC 
Mass 

Core radius 
Half-light radius 
Tidal radius 

Half- mass relaxation time (Rh) 

Binary fraction 

[Fe/H] 

Age 

A v 



a 1.72kpc 
5.9kpc 
a 63 OOOM 
0.53 pc 
2.3 pc 
21 pc 
660 Myr 

n-15% 

-1.2 

6 12Gyr 
"1.33 



R eference s: All d ata are fro m the current version of the catalogue 
of iHarrid lll996l) . except a IRicher et all (20.04) (though this is 
not always the orig inal reference for the quoted number) and b 
Hansen et all <2004h . 



The line-of-sight velocity dispersion is computed in a 
similar way. For each particle we calculate the mean square 
line-of-sight velocity (because the orientation of the trans- 
verse component of the velocity is random), and then sum 
over all particles. In this sum, each particle is weighted by a 
geometrical factor proportional to the surface density of the 
particle's shell along the line of sight. The result is an esti- 
mate of the velocity dispersion that is weighted by neither 
mass nor brightness, but only number. To check the effects 
of this, we have sometimes calculated velocity dispersions 
with various cutoffs in the magnitudes of the stars included. 
The effects are usually quite small. 

Computation of the luminosity function is the least 
problematic. We count stars in each bin in ^-magnitude, 
but lying above a line in the colour-magnitude diagram just 
below the main sequence. Again the contribution of each 
star is weighted by the same geometrical factor, and the V 
magnitude is corrected for extinction. 



2.4 Finding initial conditions 

Here we summarise our exper ience in approaching th i s prob - 
lem. A number of stu dies (e.g. lBaumgardt fc Making |2003), 
lLamers et al. I rt2005r i1 give simple formulae for the evolution 
with time of the bound mass of a rich star cluster. It is 
possible to derive similar simple formulae for the evolution 
of the half-mass radius and other quantities. There are sev- 
eral problems with inverting these formulae, however, i.e. 
using them to infer the initial parameters of a cluster from 
its present mass and radius. First, these present-day global 
parameters are quite uncertain, even to within a factor of 
two. Second, these formulae depend on the galactic orbit and 
other parameters which are equally subject to uncertainty. 
Therefore we have adopted the more straightforward but 
more laborious approach of iterating on the initial param- 
eters of our models (Tab [T] lower half); that is, we select 
values for the five stated parameters, run the model, find 
where the match with the observations is poor, adjust the 
parameters, and repeat cyclically. 

We have employed two methods to facilitate this pro- 
cess to some extent. First, we have often carried out mini- 
surveys, i.e. small, coarse grid-searches around a given start- 
ing model, to find out how changes in individual parameters 
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affect the results. Second, we have used scaling to acceler- 
ate the process, and we now describe this method in a little 
detail. 

Suppose we wish to represent a star cluster which has 
mass M and radius R with a model representing a cluster 
with a (usually smaller) mass M* and radius R*. Since two- 
body relaxation dominates much of the dynamics, we insist 
that the two clusters have the same relaxation time, and so 



R* 



M \ 1/3 Aog 7 iV 



\M*J 



log 7iV 



2/3 



where N, N* are the corresponding particle numbers. Thus 
the model of lower mass has larger radius. The tidal radius is 
scaled in the same way. Then the observational results (sur- 
face brightness profile, etc) can be computed for the model of 
lower mass and then rescaled (by appropriate factors of the 
mass and radius) to give a result for the more massive clus- 
ter (assuming that the evolution is dominated by the pro- 
cesses of relaxation, stellar evolution and tidal stripping). 
In fact we have found that this is very successful, in the 
sense that the inferred best values for the initial conditions 
change little when runs are carried out with the "correct" 
(i.e. unsealed) initial mass, certainly when the proportion of 
binaries is 10% or less. Some aspects of the evolution are 
not well described by this scaling technique. For example, 
we do not change the distribution of the semi-major axes of 
the binaries. In this way the internal evolution of the bina- 
ries is correctly modelled (provided that the binaries remain 
isolated dynamically), though their dynamical interactions 
with the rest of the system are not. In principle one could 
scale the semi-major axes in the same way as R, but then 
the internal evolution of the binaries would be altered. 



3 MODELS OF M4 

3.1 Finding initial parameters 

Our starting point was our work on the old open cluster 
M67 (G iersz fc Heggic 2008) , but with a larger initial mass 
and radius. ( Baumgardt fc Makinol (|2003T ). for example, sug- 
gested that the initial mass of M4 (NGC 6121) was of order 
7.5 x 1O 5 M0, though we used somewhat smaller values.) To 
begin with, our choice of initial tidal radius was inferred 
from the initial mass and the present-day estimates of mass 
and tidal radius given in Table [2] assuming that rt oc M 1 ^ 3 . 
At first we adopted similar values for the "concentration" 
(ft/rh) and binary fraction as in the modelling of M67, but 
found that the surface brightness profile fitted poorly (with 
too large a core) unless the binary fraction was much smaller 
(5 to 10%) and the concentration much higher. 

Before describing our best models, it is worth briefly 
mentioning one which provided a satisfactory fit to the sur- 
face brightness profiles. The fit to the velocity dispersion 
profile was tolerable, but indicated a model that was too 
massive by a factor of about 1.4. Its main flaw, however, was 
in the luminosity function, which was generally too large by 
a factor in the range 2-3. There are two reasons why this 
is interesting. One is that, for a long time, models of star 
clusters were constructed entirely on the basis of the surface 
brightness and velocity dispersion profiles. It should be re- 
alised that such models may be misleading in other ways. 



Second, this experience underscores the importance of prop- 
erly normalised luminosity functions. In other words, it is 
important to know the area of the field where the stars have 
been counted, or some equivalent representation of properly 
normalised data. Very often, the emphasis is solely on the 
shape of the luminosity function, but we have found that 
the absolute normalisation is an essential constraint. 

iBaumgardt fc Makinol (|2003T ) show that the lower mass 
function becomes flatter as the fraction of mass lost by the 
cluster increases (i.e. towards the end of its life). There- 
fore we could perhaps have improved the fit with the lumi- 
nosity function by starting with a more massive model and 
somehow ensuring a larger escape rate so as to leave a sim- 
ilar mass at the present day. Instead, we elected to change 
the slope of the low-mas s IMF from the canonical value of 
q = 1.3 (|Kroupal [2007al 'l to a = 0.9. (There is some jus- 
tification for a lower value for low-metalli city populations, 
though it has been arg ued (|Kroupall2007bl ) that there is no 
pristine low-metallicity population where the IMF can be 
inferred securely.) 



3.2 A Monte Carlo model of M4 

By some experimentation we arrived at a model which gave 
a fair fit to all three kinds of observational data; see Ta- 
ble 3 (where we co mpare with a King model developed by 
I Richer et alj (|2004f ). and Figs 1-4. It is worth noting that no 
arbitrary normalisation has been applied in these compar- 
isons between our model and the observations. The surface 
brightness profile, for example, is computed directly from 
the V-magnitudes of the stars in the Monte Carlo simula- 
tion, as described in Sec l2.3l In the construction of a King 
model, by contrast, it is often assumed that the mass-to-light 
ratio is arbitrary. 



3.2.1 Surface brightness 

While the overall surface brightness profile is slightly faint, 
the most noticeable feature of FigfJ] is that the model has 
a somewhat smaller limiting radius than the observational 
data. The reason for this is explained in iGiersz fc"~f lcggic 
(2008): in short, we impose a smaller tidal radius than the 
nominal tidal radius, in an iV-dependent way, which is in- 
tended to ensure that the overall rate of escape from the 
model behaves in the same way as in an JV-body simulation. 
Another point to notice is the disagreement between the 
total luminosity of our model and tha t of the King mode l 
quoted in the final column of Table [3] iTrager et al.l (|l995T ) 
give an analytic fit to the surface brightness profile, and we 
have checked that the integrated value is close to ours. 

The data for the Monte Carlo model in Table [3] would 
be consistent with a galactocentric radius of about 1.7kpc, 
in an isothermal galaxy model with a circular velocity of 
220km/s. While this is certainly much smaller than its cur- 
rent galactocentric distance, a small value was also found 
(using a similar arg ument and publ i shed v alues of the mass 
an d tidal radius) by van den Berghl (|l995h . The orbit given 
bv lDinescu et all (|l999T ) has a still smaller perigalactic dis- 
tance, the galactocentric distance varying between extremes 
of 0.6 and 5.9 kpc. 
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Table 3. Monte Carlo and King models for M4 



Quantity 


MC model 
(t = 0) 


MC model 
(t = 12Gyr) 


King model 
(Richer et al. 2004) 


Mass (Mq) 


3.40 x 10 5 


4.61 x 10 4 




Luminosity (Lq) 


6.1 x 10 6 


2.55 x 10 4 


6.25 x 10 4 


Binary fraction f b 


0.07 


0.057 





Low-mass MF slope a 


0.9 


0.03 


0.1 


Mass of white dwarfs (Mq) 





1.81 x 10 4 


3.25 x 10 4 * 


Mass of neutron stars (Mq) 





3.24 x 10 3 




Tidal radius rt (pc) 


35.0 


18.0 




Half-mass radius (pc) 


0.58 


2.89 





*: this is the quoted mass of "degenerates" 
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Figure 1. Surface brightness p rofile of our Monte Carlo 
model, compared with the data of lTraeer et al.l dl995h . 
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Figure 2. Velocity dispersion prof ile of our Monte Carlo 
model, compared with the data of I Peterson et all (|l995). 




Figure 3. Luminosity function of our Monte C arlo model at 
the m edian radius of the innermost annulus in R icher et al.l 
(2001). compared with their data. 



Outermost annulus 




Figure 4. Luminosity function of our Monte C arlo model at 
the m edian radius of the outermost annulus in iRicher et al.l 
(2004), compared with their data. 



6 M. Giersz and D. C. Reggie 



3.2.2 Velocity dispersion profile 

This is illustrated in Fig[2l where i t is co mpared with the 
observational data of lPeterson et all (j 19951 ) . The shortfall at 
large radii is of doubtful significance. 



3.2.3 Luminosity Functions 

These are shown for our model at the medi an radius of 
the in nermost and outermost fields observed bv lRicher et al.l 
(2004). The disagreement in the inner luminosity function 
at faint magnitudes may be attributable to the fact that 
the theoretical result assumes 100% completeness, while the 
observational data are uncorrected for completeness. A plot 
of the comple t eness correction in the outer field is given by 
IHansen et al.l (|2002l . Fig. 3) and discussed above in Sec l2.3l 
The mismatch between model and data is smaller in the 
outer field, but from the discussion above it would seem 
that the mismatch in the faintest bin may be too large to 
be accounted for by the estimated value of the complete- 
ness correction. The error bars in the last two observational 
points on this plot (not shown) almost overlap, and much 
of the mismatch may be simply sampling uncertainty. It is 
worth comparing th e multi-mass King model constructed by 
iRicher et all l|2004h , which is also problematic in the outer- 
most field. 



3.2.4 Core collapse 

Fig[5] shows the evolution with time of the theoretical core 
radius. There is an early period of very rapid contraction, 
associated with mass segregation, followed by a slower reex- 
pansion, caused by the loss of mass from the evolving mas- 
sive stars which are now concentrated within the core. 

Our most surprising discovery from our model of M4 is 
the subs equent behaviour. M4 is classified as a King-profile 
cluster (|Trager et al.lll993h . and such clusters are usually 
interpreted as being clusters whose cores have not yet col- 
lapsed. But the plot of the theoretical core radius (Fig[5] re- 
veals that the model exhibited core collapse at about 8Gyr. 
Subsequently its core radius is presumably sustained by bi- 
nary burning. Even non-primordial binaries may be playing 
a role here. To the best of our knowledge it has not pre- 
viously been suggested that M4, is a post-collapse cluster, 
thoug h on statistical grounds lDe Marchi. Paresce. fc Pulond 
(|2007l ) have suggested that some King-type clusters have al- 
ready collapsed. This issue is often approached by reference 
to the half- mass relaxation time (Table [2] but for M4 this is 
not short enough to be decisive. 

The presence of radial colour gradients is correlated 
with the presence of a non-King surface brightness profile 
and hence with core collapse. We have computed the colour 
profile of our model, and find that it is nearly flat except 
for the influence of one or two very bright stars within the 
innermost few arcsec. Most values lie around 0.65 in B — V , 
which i s muc h less than the global value of 1.03 given in 
iHarrid (|l996l . January 2008). On the other hand roughly 
similar mismatches between observations are found in other 
clusters (e.g. M30, (jPiotto. King, fc Diorgovskilll988l) ). 




2000 4000 6000 8000 10000 12000 

Time (Myr) 

Figure 5. Theoretical core radius of our Monte Carlo model. 



3.2.5 The colour-magnitude diagram 

Figure [6] shows the colour- magnitude diagram of the model. 
This is of interest, not so much for comparison with ob- 
servations, but for the presence of a number of interesting 
features. The division of the lower main sequence is simply 
an artifact of the way binary masses were selected (a total 
mass above O.2M0 and a component mass above O.IMq.) Of 
particular interest are the high numbers of merger remnants 
on the lower white dwarf sequence. There are very few blue 
stragglers. Partly this is a result of the low binary frequency, 
but it is also important to note that some formation chan- 
nels are unrepresented in our models (in particular, collisions 
during triple or four-body interactions, though if a binary 
emerges from an interaction with appropriate parameters, it 
will be treated as merged.) These numbers also depend on 
the assumed initial distribution of semi-major axis, which is 
not yet well constrained by observations in globular clusters. 

Note in Fig(6] that there are some whit e dwarf-main 
sequen ce binaries (below the main sequence). IRicher et al.l 
(|2004h drew attention to the possible presence of such bina- 
ries in their colour magnitude diagrams of this cluster. 



3.2.6 Binaries 

Photometric binaries are visible in Fig[S] and t hese are com- 
ared with observations in the inner field of IRicher et al.l 
20041 1 in FiglTl In this fi gure, the model histogram has been 
normalised to the same total number of stars as the obser- 
vational one. We made no attempt to simulate photometric 
errors, but the bins around abscissa = -0.75 suggest that the 
binary fractions in the model (which is under 6% globally; 
see Table EH) and the observat ions are comparable. We note 
here that Richer et al.l |2004), using this same data, con- 
cluded that the binary frequency was approximately 2% in 
the innermost field. But they also note that the measured 
frequency of "approximately equal-mass binaries" is 2.2%, 
and we consider that the balance between this number and 
our binary fraction can be made of binaries with companions 
whose masses are more unequal. 
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Figure 6. The colour-magnitude diagram at 12Gyr. Green: single 
stars; red: binaries; blue pluses: collision or merger remnants. 



Radial velocity binaries should also be detectable in M4, 
and will be part of a separate investigation. 

Now we consider the distributions of the binaries, as 
predicted by the theoretical model. Fig(8] shows that bina- 
ries have evolved dynamically as well as through their in- 
ternal evolution. In particular the softest pairs been almost 
destroyed. 

By 12 Gyr the binaries exhibit segregation towards the 
centre of the cluster, but perhaps in more subtle ways than 
might be expected (Figs l9llOTl . When all binaries are consid- 
ered, there is little segregation relative to the other objects 
in the system. (Most binaries in our model are of low mass.) 
But if one restricts attention to bright binaries, which we 
here take to mean those with My < 7 (i.e. brighter than 
about two magnitudes below turnoff ), the segregation is very 
noticeable fFig |10|l . with a half- mass radius smaller by al- 
most a factor of 2 than for bright single stars. Still, bright 
binaries are not nearly as mass-segregated as neutron stars 
(Fig[9]), which, incidentally, receive no natal kicks in our 
model. 

The history of the binary fraction is, in effect, given 
in Fig |13l In the first Gyr this falls steadily from the initial 
value of 0.07 to about 0.06, and it remains close to this value 
until the present day. 



3.2.7 Escape velocity 

We have already referred to the fact that, in our model, 
neutron stars receive no natal kicks. Escape is of course gov- 
erned by the escape speed, and this is of much interest in 
connection with the possibility of retaining gas from the first 
generation of rapidly evolving stars. For this reason we plot 
the central escape velocity in Fig. 1111 The remarkably high 
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Figure 7. Histogram of V-oSset from the main sequence, 
compared with the corresponding data from the innermost 
annulus studied by Richer et al. (2004) See text for details. 
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Figure 8. Distribution function of the semi-major axes of the 
binaries at OGyr and 12Gyr. Units: solar radii. 



value in the first few tens of millions of years, if valid for 
M4, could have interesting consequences for the early evo- 
lution of the cluster and its stars. It draws attention to the 
very high initial density of our model (Table whose av- 
erage value within the half-mass radius is 2 x lO 5 M0/pc 3 . 
The central density is about lO 6 M0/pc 3 , about an order of 
magnitude larger tha n in the central yo ung cluster in the 
HII region NGC 3603 (|Stolte et alj|2006h . In the absence of 
local young star clusters as massive as our M4 progenitor, it 
is difficult to be sure whether our initial model is implausibly 
dense. 



3.2.8 Degenerate components 

We have already mentioned the spatial segregation of the 
population of neutron stars, and the presence of a number 
of degenerate binaries at the present day. Now we consider 
the historical evolution of this population over the lifetime 
of the cluster so far, according to our model. 

Fig |12l shows the evolution o f these populations. It has 
already become well established (|Vesperini fc Heggielll997l ; 
iHurlev fc Sharall2003D that white dwarfs account for an in- 
creasing proportion of the mass of globular clusters, and in- 
deed it is of order 39% at the present day, according to this 
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Figure 9. Radial distribution functions at 12Gyr. Units: par- 
sees. The key identifies the class of object included. The dis- 
tributions of white dwarfs and binaries are almost identical. 
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Figure 12. Mass of white dwarfs, and of neutron stars and black 
holes together. 
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Figure 10. Radial distribution functions at 12Gyr, showing 
the extent of segregation between bright single and bright 
binary stars (My < 7). Also shown for comparison is the 
distribution for all objects, as in Fig[9] Units: parsecs. 
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Figure 13. Numbers of all stars, white dwarfs, neutron stars and 
black holes. 
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Figure 11. Central escape velocity as a function of time. 



model. The proportion of neutron stars is almost certainly 
excessive, because of our assumption of complete retention. 

In order to separate two of the components in this fig- 
ure (black holes and neutron stars), we show in Fig. 1131 the 
numbers of these degenerate stars, along with the numbers 
of all stars and all binaries. This shows that no stellar-mass 
black holes are expected to be present in M4 now. 



4 CONCLUSIONS AND DISCUSSION 

In this paper we have presented a Monte Carlo model for 
the nearby globular cluster M4. This model includes the 
effects of two-body relaxation, evaporation across the tidal 
boundary, dynamical interactions involving primordial and 
three-body binaries, and the internal evolution of single stars 
and binaries. By adjustment of the initial parameters (total 
mass, tidal radius, initial mass function, binary fraction) we 
have found a model which, after 12Gyr of evolution, leads to 
a model with a surface brightness profile, velocity dispersion 
profile, and luminosity functions at two radii, all of which 
are in tolerable agreement with observational data. It also 
leads to a number of photometric binaries roughly consistent 
with observation. 
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This model has a current mass of 4.6 x 1O 4 M0, a M/Lv 
ratio of 1.8, a binary fraction of almost 6%, and an almost 
flat lower mass function (Table [3]). Almost 40% of the mass 
is in white dwarfs, and about 7% in neutron stars, though 
our model assumes a 100% retention rate. They are strongly 
mass-segregated to the centre of the cluster. There are no 
stellar-mass black holes left in the cluster. The binaries have 
experience significant dynamical evolution, almost all the 
soft pairs having been destroyed. The binary population as 
a whole is only slightly segregated towards the centre, but 
there is more evident segregation of bright binaries (such 
as those that would be more readily observed in a radial 
velocity search.) 

The most significant new result from our model is the 
implication that M4 is a core collapse cluster, despite the 
uncollapsed appearance of the surface brightness profile. 

Now we discuss a number of shortcomings of our model 
and other issues related to this study. 

(i) Uniqueness: Though we have arrived at a broadly sat- 
isfactory model, it is not at all clear how unique it is. Cer- 
tainly we were unable to construct models with much larger 
numbers of primordial binaries, or with significantly larger 
initial radius, as such models produced an insufficiently con- 
centrated surface brightness profile. Furthermore, the fact 
that the lower slope of the initial mass function is less steep 
than the canonical value of 1.3 should not be taken to imply 
that we have inferred the initial value. 

(ii) Fluctuations: we have noticed that runs differing only 
in the initial seed of the random number generator can give 
surprisingly different surface brightness pr ofiles. In broa d 
terms this confirms the important finding of lllurlevl (|2007h . 
though we have not yet established (as he did) that the 
presence or absence of black hole binaries is the underly- 
ing mechanism. We shall return to this issue at appropriate 
length in the next paper in this series, on the cluster NGC 
6397. 

(iii) The Initial Model: The initial model is astonishingly 
dense, the central density being of order lO 6 A/0pc -3 . It 
also has to be realised that we are imposing initial con- 
ditions at a time when the residual gas from the birth 
of the cluster has already dispersed. Various properties of 
a star cluster, including its binary population and mass 
function, may change significantly during the phase of 
gas expulsion, which further undermines the power of our 
model to establish the initial conditions. Finally, we have 
assumed no initial mass segregation, even though recent 
work suggests that this should be present already before 
the cluster has assembled into a roughly sph erical object 
^McMillan. Vesperini. fc Portegies ZwartlfeoOTl ) . This is an 
aspect of the modelling that could be readily improved. 

(iv) The early evolution of the model: It is worth drawing 
attention here to the high initial escape velocity from the 
centre of our model, in the first few tens of millions of years. 
Clearly this is dependent on the small initial radius. The 
initial high density is also responsible for the very rapid 
initial evolution of the core. 

(v) Imperfections of the modelling: Several important im- 
provements need to be made to our technique. 

(a) At present few-body interactions are handled with 
cross sections. The problem is not simply that these are 
not well known, especially in the case of unequal masses; 



it also means that we are unable to determine if a collision 
occurs during a long-lived interaction. For this reason, we 
have said almost nothing about collision products (e.g. 
blue stragglers) in this paper. This limitation could be 
overcom e by direct integr a tion o f the interactions, as is 
done by iFregeau fc Rasiol f2007|) in their version of the 
Monte Carlo scheme. 

(b) Long-lived triples are neglected in the model at 
present. Th ese are commo nly produced in binary-binary 
encounters l|Mikkolalll98i ). and it is desirable to include 
these as a third species (beyond single and binary stars). 
Their observable effects may be small, but of course there 
is one intrigui ng example in the ve ry cluster we have fo- 
cused on here (|Thorsett et al.lll993h . 

(c) We assume that the tide is modelled as a tidal cut- 
off. We have taken some care to ensure that the tidal 
radius is adjusted so that the model loses mass through 
escape at the same rate as an TV-body model would, 
if immersed in a tidal held with the same tidal radius. 
This means, however, that the effective tidal radius of the 
model is somewhat too small. A better treatment of a 
steady tide may be possible, without this drawback. 

(d) We assume that the tide is steady, something which 
is not true for M4. The effects of tid al shocks have been 
studi ed by a number of authors (e.g. iKundic fc Ostrikerl 
1995), and it would be possible to add the effects as 
another process altering the energies and angular mo- 
menta of the stars in the simul ation. On the other hand 
lAguilar. Hut, fc Ostrikerl (|l988T l found, on the basis of a 
simple model, that tidal evaporation was the dominant 
mechanism in the evolution of the mass of M4 at the 
present day, and that other factors (such as disk and bulge 
shocking) contributed at a level les s than 1%. 

(e) Rotation: it has been shown (|Kim et al.ll2004l ) that, 
to the extent that rotating and non-rotating models can 
be compared, rotation somewhat accelerates the rate of 
core collapse. Rotation is hard to implement in this Monte 
Carlo model, however. 

(f) The search for initial conditions is still very labo- 
rious, and we are constantly seeking ways of expediting 
this. Our most fruitful technique at present is the use of 
small-scale models which relax at the same rate as a full- 
sized model. Not all aspects of the dynamics are faithfully 
rendered by a scaled model, but the technique appears to 
be successful as long as the fraction of primordial binaries 
is not too large. 

Despite these caveats and shortcomings, it is clear that 
the Monte Carlo code we have been developing is now an 
extraordinarily useful tool for assessing the dynamics of rich 
star clusters. For the cluster M4 it has led to the conclusion 
that M4 can be explained as a post-collapse clusters. It also 
provides a wealth of information on the distribution of the 
binary population, which will be important for the planning 
and interpretation of searches for radial velocity binaries, 
now under way. 

A-body models will eventually supplant the Monte 
Carlo technique, but at present are incapable of providing a 
star-by-star model of even such a small cluster as M4. They 
can of course provide a great deal of general guidance on 
the dynamical evolution of rich star clusters, and they un- 
derpin the Monte Carlo by providing benchmarks for small 
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models. Even when iV-body simulations eventually become 
possible, Monte Carlo models will remain as a quicker way 
of exploring the main issues, just as King models have con- 
tinued to dominate the field of star cluster modelling even 
when more advanced methods (e.g. Fokker-Planck models) 
have become available. 
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